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ABSTRACT 


Dual 2ind multi-spool rotors are being increasingly used in large output compressors and 
turbines. In such applications the rotor size becomes large, which leads to increase in 
shaft flexibility and associated problems. In order to overcome this constraint and to get 
additional advantages single rotors are converted to twin or multi-spool rotor systems, 
whereby shafts are made hollow to support other shafts through inter shaft bearing(s) The 
focus of the present study is to model such shaft-disc-bearing configurations for 
rotordynamic analysis. Transfer matrix approach has been adopted to develop the 
analytical model and programming is done using MATLAB. In addition Finite Element 
models have been developed using NISA package. Free vibration analysis, including 
estimation of natural frequencies, mode shapes, critical speeds for backward and forward 
whirls and unbalance response analysis including whirl orbit computation and beat 
phenomenon have been carried out. Standard problems available in literature are 
considered and results obtained through both procedures are compared to validate the 
models. Transfer matrix and FEM models have also been developed for a typical two- 
spool multi-stage aeroengine rotor using sample geometric, material and bearing data. 
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CHAPTER 1 


INTRODUCTION 


More than ninety percent of machinery used in industry can be categorized as rotating 
machines. Some of them are steam and gas turbines, compressors, turbogenerators, 
aeroengines etc. These are typical multimass rotor installations and comprise of 
mechanical elements like shafts, discs, blades, seals, couplings, gears etc. The basic 
model is constituted of a shaft carrying one or more than one disc. The arrangement may 
be supported in rolling element or in flxxid film bearings. 

In case of large output compressors and turbines, the rotor size becomes large, which 
leads to increase in shaft flexibility and other associated problems. In order to overcome 
this constraint and to get additional advantages a single rotor is converted to a twin or 
multi-spool rotor system, whereby shafts are made hollow to support other shafts, 
through inter shaft bearing(s). Such rotor systems are widely used in aeroengines, where 
various compressor arid turbine stages (LP, IP, HP) are combined together to form a set 
of inner and outer rotors. A dual compressor aeroengine utilizes two separate 
compressors, each with its own driving turbine. One of the principle advant^es of 
split compressor in aeroengine is greater flexibility of operation without danger of stall. 
The low-pressure compressor can oj^rate at the best speed for accommodation of the low 
pressure, low temperature at the forward part of engine. The high-pressure compressor is 
speed governed to operate at the proper speed for most efficient performance. Use of dual 
compressor makes it possible to attain pressure ratio of more than 13:1, where single 
axial compressor produce pressure ratio of 6 or 7:1 . 

This type of rotor system minimiz es shaft dejections caused by rotor unbalance and 
improves engine efficiency, p^ormance ami peliabilky by eliminating the static support 
structure in the aerodynamic flow. It also fadlitafe^ easy mounting with engme casing 
and provides compactness to the overall ermtoe. 



However, the intershaft bearing(s) becomes a source of cross-excitation between the two 
shafts. This cross-excitation makes the mathematical treatment of dual rotor system 
different from that of single rotor. Cross-exciting vibration between inner and outer shafts 
is effected through the intershaft bearing and thus the dynamic response of one rotor also 
depends upon the dynamic behavior of the other. 

The present study is an attempt in rotor dynamic modeling and analysis of two/mxilti- 
spool rotors. Rotor dynamics has been an area of extensive investigation in the past. Two 
distinct approaches for analysis of multimass rotors have been the continuous systems 
approach and the discrete systems approach. In continuous systems approach the mass 
and stiffness properties are taken to be distributed continuously along the shaft. This 
method is used in finite element method with component mode s 5 mthesis to solve dual 
rotor or large system problem. In discrete systems approach mass and stiffiiess properties 
are discretized and matrix methods are generally employed to obtain system response. 
The finite element method has a distinct advantage in assembling the elemental equation 
with out any recourse to intershaft bearing conditions as in the case of transfer matrix 
method. Dimentberg (1961) gave an exhaustive treatment to the transverse vibration of 
rotating shafts through a continuous system approach. Tondl (1965) analyzed various 
aspects of rotor dynamics, like critical speeds, stability, self excited vibrations, fluid film 
bearing characteristics etc. Another text on the subject is by Rao (1996) who bases his 
treatment on the discrete system approach, specifically transfer matrix methods of 
Myklestad (1944) and Prohl (1945). Some salient studies on various aspects of rotor 
dynamics are - Jeffoct on unbalance response of single mass rotors, Kikuchi (1970), on 
vibration of multi-disc and multi-bearing rotor systems; Green (1958) and Carnegie 
(1964), on gyroscopic effects; Smith (1969), Limd (1965) and Morton, Johnson and Wale 
(1988) on fluid film bearings. Goodydn (1955) has carried out extensive experimental 
studies for the determination of bearing coefficients, which are eight in number. He has 
presented his results in the form of grjq)hs and diarts. Lalanne (1985) has reviewed 
vibration problems of jet engines. Childs (1976) developed a modal simulation model 
based on eigendata at zero running ^)eed, accomfting for gyroscopic effects, bearing 
damping and nonlinearities, structural modal damping and concentrated damping due to 
squeeze film dampers. Glasgow and Nelson (1979) ^plied component mode synthesis 
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and showed that significant reduction in the size of the problem is achieved while still 
retaining the essential dynamic characteristics of the lower modes. Other analysis using 
finite element method have been done by Berthier et al. (1985) and using component 
mode synthesis by Li and Gunter (1982). Zeng and Hu (1988) use gyroscopic mode 
synthesis technique for multishafl rotor-bearing case system. Squeeze film dampers play 
an important role in multi spool aeroengine rotor system. Studies on stability aspects 
using intershaft squeeze film damper in two spool rotor system have been carried out by 
several authors, Hibner, Bansal and Buono (1978), Li and Hamilton (1986) and by El 
Shafei (1991). Hibner et al. (1977, 1978) investigated the stability behavior of dual rotor 
system with squeeze film damper and suggested means for controlling the instability by 
changing the stiffness of the intershaft bearings by adding the spring in parallel. The 
transfer matrix method has been applied in various forms to solve multi spool rotor 
system problems. Hibner (1974) applied a unique transfer matrix method to an idealized 
equivalent engine system for predicting vibratory response accounting for nonlinear 
viscous damping effects. Another analysis has been by Bansal and Kirk (1975) for 
stability and damped critical speeds. Kazao and Gunter (1989a) have also studied energy 
distributions in the rotor system and have shown its utility in designing effective bearings 
and dampers. Gupta et al. (1988) presented a theoretical formulation to determine the 
steady state unbalance response of dual rotor system with a flexible intershaft bearing 
using extended transfer matrix method, where the extended transfer matrix assumed a 
dimension of 33x33. Transfer matrix impedance <x)upling method developed by Taiping 
(1988) is useful to find eigensolutions of multi-spool rotor system. 


Present Work 

The focus of the present work is to develop models for dynamic analysis of two and 
multi-spool rotors. The modeling involves (i) free vibration analysis - including 
estimation of natural firequencies, mode shapes, critical speeds for backward and forward 
whirls and (ii) unbalance response analysis including whirl orbit computation and beat 
phenomenon. Models have been developed ftnough (i) transfer matrix approach with 
programming being done in MATLAB and (ii) Finite Element procedure using NISA 
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FEM package version 7.0. Standard problems available in literature are considered and 
results obtained through both procedures aire compared to validate the models. Features 
such as the nature of the steady state orbit, cross excitation between the rotors, critical 
speeds and mode shapes are discussed. Free vibration modeling procedure is discussed in 
Chapter 2, while unbalance response analysis is discussed in Chapter 3. In addition the 
modeling exercise has been extended to analyze a typical two-spool multi-stage 
aeroengine rotor. Rotordynamic model is developed from sample geometric, material and 
bearing data. This is described in Chapter 4. Detailed parametric studies could not be 
carried out. Preliminary results on the free and residual unbalance response are presented. 




CHAPTER 2 


CRITICAL SPEEDS OF DUAL ROTORS 


The first step in the dynamic analysis of any rotor-bearing system comprises of 
computation of critical speeds and mode shapes. Unbalance response analysis and 
consequent dynamic stress analysis is carried out after necessary critical speed 
calculation and free vibration analysis has been accomplished. In the present study, 
dual rotors are modeled through two separate procedures, namely - (i) transfer matrix 
method of Myklestad and Prohl (1944) and (ii) finite element technique, using 
standard available software. The modeling procedures and results are described in this 
chapter. 

2.1 Transfer Matrix Method 

The transfer matrix method, for far-coupled systems, suggested by Myklestad and 
Prohl has proved to be an effective tool for modeling rotor-bearing systems for 
bending vibration analysis. State vectors defined at various discrete stations along the 
rotor length are related through Field and Point matrices and boundary conditions are 
applied to obtain the frequency equation, which can be solved through any numerical 
technique. 


2.1.1 Field Matrix 

An n mass rotor system is shown in Fig 2.1(a). The rotor discs are located at stations 

1,2, n. The rotor is supported at stations 0 and n+1. The masses are taken to 

be lumped with gyroscopic inertia neglwjted. Defining the state vector at station i as 

f-wl 


e 



> 


( 2 . 1 ) 



where w is deflection; 0 is the slope; denotes the bending moment and F_ 


denotes shear force. The state vector can be defined to the left or rights of station i, as 
shown in Fig 2.1(b). /,■ is the shaft element length and w, is the mass located at the zth 
station in Fig 2.1(b). The relations across the ith field can be written as firom Fig 2.2 
(Rao, 1996) 
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Defining a Field Matrix 
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equation (2.3) can be symbolically written as 

{5}f =[Fl{5}f_i 


(2.2) 


(2.3) 


(2.4) 


(2.5) 


2.1.2 Point Matrix 

The force and moment balance relations across the itii mass (Fig 2.3) provide the 


following matrix equation 
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( 2 . 6 ) 


Defining the /th Point Matrix as 


[n = 


1 

0 

0 

mp^ 


0 0 0 
1 0 0 
0 1 0 
0 0 1 


(2.7) 


equation (2.6) can be symbolically expressed as 

{5}f= [/>],. {S>f (2.8) 

2.1.3 Overall Transfer Matrix and Frequency Equation 

Employing equations (2.5) and (2.8), the matrix relationship across stations 0 and n+1 
can be written as 

Wf =mi{s>o 

{S}f =[f]i[F]i{S}o 

-[FbWf =[F]2[/>ht^')l{-S'}0 (2-9) 

=[F)„+,[P]„[F]„[?]„., [F],{S}0 

Defining the product of all field and point matrices, in the order given in the above 

C' ' ' 

equation, as the overall transfer matrix [(/) one gets 

{S>™+l=[C^]{S}o , (2-10) 



Since field and point matrices are both of order 4. the overall transfer matrix is also 
4x4 in size and equation (2. 1 0) can be written in expanded form as 
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( 2 . 11 ) 
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Equation (2.11) is employed, after substitution of the boundary conditions 
appropriately, to obtain the firequency equation of the rotor and consequent 
computation of its critical speeds. For example, for a rotor supported in rigid bearings 
and modeled as a simply supported system, the boundary conditions are w= 0 and 
My = 0, at stations 0 and n+\, which upon substitution in equation (2.1 1) yield the 
following frequency equation 


“12 “14 

|“32 “34 


:0 


( 2 . 12 ) 


An appropriate root searching technique can be employed to search for the values of 
critical speeds pi in equation (2.12). Once the critical speed is determined, the 
corresponding mode shape is determined by assuming a unit slope at station 0, i.e., 
^ 0=1 which gives Fj o = -Wi 2 /“i 4 thtis the state vector at station 0 is 

completely defined. Employing the relations (2.5), (2.8) one can proceed to determine 
the state vectors at all stations and get the mode shape corresponding to the critical 
speed determined. 


2.1.4 Gyroscopic Effects 

Gyroscopic couple changes the bending moment equation across the mass, which is 
not considered earlier. The disks on the rotor contribute to the gyroscopic couple. For 
discs located at nodal points, the gyroscopic effect can be significant, especially in the 
case of overhang rotors, with a disk at the free end. Considering the spinning of the 
rth disk of the rotor, if the angular velocity is co and the transverse and polar moment 
of inertia of the disk are Ij- and Ip respectively, the precessional motion ^and in 
the x-z and x-y planes respectively, are governed by (refer Fig.2.4 a, b) 

.11 ■ 





(2.13) 


^yi ~ ^ yi I pO)^i+ Ij 
M^i = M^i - 1 pO}0i+ Ij (t)^ 

which modify the point matrix to 

1 0 0 0 

0 1 0 0 

0 {Ip p)p 1 0 

THiP^ 0 0 0 

It is to be noted that forward synchronous whirl of the rotor would correspond to the 
case G) = p, while backward synchronous whirl corresponds to co = -p. The effect 
of gyroscopic couple is to stiffen the rotor and raise the critical speed. Due to 
gyroscopic effect the natural frequency of rotor system changes with rotating speed. 

2.1.5 Intershaft Coupling 

The traivsfer matrix procedure employing point and field matrices described above 
with appropriate modifications for gyroscopic effects is an efficient tool for finding 
the natural frequencies of single-shaft / spool rotor system. However, in case of multi- 
shaft / spool rotor system the dynamic behavior of one rotor is affected by another 
due to the presence of intershaft bearings and this makes its mathematical analysis 
different 

Figure 2.5 shows the scheme of a typical multi-spool rotor system with several 
intershaft bearings. It consists of N rotors and M intershaft bearings. The system is 
divided into N+M subsystems - N rotor systems arKi M inter shaft bearings 
subsystems, i.e., rotor 1, rotor 2, ...rotor N and inter-shaft bearings with impedances 
Zbj, Zb 2 :■ Zbm- On separation of the connections of every intershaft bearing, the 
constraint condition will be a p^ of external force vectors Rjj, added on the 
connecting point j, in the y and z directions respectively. The transfer matrix method 
is now employed for ftie rotor • sub-systdns. For subsystem rotor i, the boundary 
condition will be 

vl3 


(2.14) 




(2.15) 


m;, = M'o = 0, 

My}^ = = 0, 

K,=n=o 

The relationship between the state variables on section 0 and the left side of station 1 
in zth rotor will be 

{«};''= [f/liWi (2,17) 

The state vector on the right side of station 1 will be 

{S}[^ = {- w, e, My,V^+R^, V, (!>, -Vy + Ry (2. 1 8) 

The relationship between the state vector on the sections of the right side of station 1 
and left side of station 2 will be 

(2.19) 


Rotor 1 



Figure 2S General multi-spool rotor system scheme 


Similarly, the state vector on the section of liie ri^t side of station 2 will be: 
{S)f ={-w,0,M,,r, +R„v4,M,-V, 


14 


( 2 . 20 ) 



The relationship between the state vector on section of the left and right sides of the 
segments, s-k, will be: 


Wl =[UMS}f 

where 

{S}f = {-w,0,My,V,+R„v,^,M,,-Vy+R^}J 


( 2 . 21 ) 


( 2 . 22 ) 


According to the boundary conditions 
obtained: 

Kk = iuUxj){S}UJ)=o 

y=i 




(2.15, 2.16), the following equations can be 


(2.23) 


K, = l:c/'(7,y){5}'^(y)-o 

M 

= l:t^l(8,y){5}'^(;) = o 

M 


(2.24) 


Substituting equations (2.17)-(2.22) in to equation (2.23) and (2.24), the following 
four homogeneous ordinary equations will be obtained 

A(4i - 3,4/ - 3)wi + A(4i - 3,4/ - 2)0^ + A(4i - 3,4/ - l)v^ + A(4i - 3,40^^ 

+ A(4i-3,4^ + l)Jiyi + + ^(4/-3,4iV + 5).K,^ + J(4/-3,4iV + 5 + l)i?,i 

+ + A(4i - 3,4N + 2s) R^ =0 

Ai4i - 2,4/ - 3)w^ + A{4i - 2,4/ - 2)^^ + A(4/ - 2,4/ - l)vj +A(4i - 2,4i)^o 

+ A(4i - 2,4N + l)Ryi + + A(4/ - 2,4iV + + A(4i - 2,4iV + ^ + l)R,i 

+ A(4i — 2,4N •^'2s)R^g —0 

.4(4/ - 1,4/ - 3)wo + A(4/ - 1,4/ - 2>9o + ^<4/ - 1,4/ - l)vo + ^(4/ - l,4/)^o 

+ A{4i - 1,4 A" + \)Ryi + + A(4i- l,4iV + s)Ry^+ A(4i -l,4N + s + l)R^i 

+ + A(4/ - l,4iV + 25)i?^ = 0 
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AiAi,Ai - 3)w^ + A{Ai,Ai - 2)^^ + A{Ai,Ai - l)i/^ + A(Ai,Ai)(l>l^ + AiAi,AN + + 

+ A{Ai,AN + s)Ry^ + A{Ai,AN + 5 + + + A{Ai - 3,4iV + 2s)R^ = 0 

(2.25) 


where A(i,j) is a function of the concerned transfer matrices [t/]l,[t /]2 [f^]* • For 

N rotor subsystem, 4A^ equations are obtained. There are, however, AN+2M unknown 
state variables wl,el,v\,(l>l, ,9^ ,v^ ,(l>^ ,Ry^,...RyM,R,i R,m- The 2M 

supplemental equations for searching for all these variables can be obtained by 
impedance matching at the coimecting points. Considering an intershaft bearing 
connecting two-rotor subsystem, as shown in Fig 2.6, the following equation must be 
satisfied 







(2.26) 


Rotor i 


Rotor i+1 



Figure 2.6 Bearing Impedance 

Deflections , in equation (2.26), can be expressed as the functions of 

Wq,9q, and the concerned transfer matrices [17]| [i7]y and 

so tiiat, equation (2.26) can be written as the following set of 

equations 

A(AN + j,Ai - 3)w^ + A(AN + j,Ai - 2)91 + - l)v^ + A(AN + j,Ai)(j>'Q 

+ A{AN + 7,4/ + l)w^'"‘ + A{AN + jM + 2)9^^'^ + 4(4A^ + j,Ai + 3)vi"'‘ + 

A{AN + j,Ai + A)4^^^+A{AN + j, AN + \)R^y ^A{AN + j, AN + j)R^ 

+ A{AN + j,AN-^M + \)Ryx+ + A{AN + j,AN + M + j)Ryi 


(2.27) 



A{AN + M + j,Ai - 3)w^ + A{AN + M + j,Ai - 2)/9i + A{AN + Af + j,Ai - 1 )v'q + 

AiAN + M + J,Ai)4>o + + M + j,Ai + + A(AN + M + j,Ai + + 

A(AN + M + j,Ai + 3)vj+' + Ai_AN + M + j,Ai + A)(^l*' + AiAN + M + j,AN + 1)/?^, 

+ + A(AN + M + j,AN + j)R^j + AiAN + M + J,AN + M + + + 

AiAN + iW + j,AN + M + j)Ryj = 0 

(2.28) 

For M intershaft bearings, 2M homogeneous equations can be obtained. Now a set of 
homogeneous ordinary equations will consist of AN+2M equations according to 
equations (2.25), (2.27) and (2.28) and can be represented in matrix form as 
[^]{5} = 0 (2.29) 

where 

{5} = {wJ,0i,V^,4 RyMf (2-30) 

Equation (2.29) can be solved through an appropriate routine to obtain eigenvalues 
and eigenvectors. 

2.1.5 Case Study 1 

The application of the transfer matrix method, described above, is illustrated through 
the following example. An undamped dual rotor containing two isotropic intershaft 
bearings, as shown in Fig.2.7 is considered. Each shaft carries two discs. An exact 
solution of this problem has been obtained by Taipeng (1988). He also employed the 
method of component mode sjmftiesis and compared the results with the exact 
solution. 






The division into subsystems - rotor 1, rotor 2 and intershaft bearings kg and kn is 
shown in Fig.2.8 



Figure-2.8 Sub-systems of the dual rotor - Case Study 1 


Table 2.1 Dual rotor system data -Case Study 1 


Station 

No. 

Mass 

(Kg) 

2 

(Kg-m") 

K 

(N/m) 

/ 

(m) 

1 

00.325 

0.0 

2.5x10^ 

0.075 

2 

10.450 

0.113 

- 

0.075 

3 

00.935 

0.0 

- 

0.225 

4 

01.321 

0.0 

- 

0.225 

5 

00.925 

0.0 

- 

0.075 

6 

08.450 

0.090 

- 

0.075 

7 

00.305 

0.0 

1.5x10’ 


8 

00.193 

0.0 

1.2x10’ 


9 

07.780 

0.071 

- 

0.075 

10 

00.373 

0.0 


0.150 

11 

06.280 

0.043 

- 

0.150 

12 ' 

00.173 

0.0 

l.OxlO’ 

0.075 


The rotational speeds of rotor 1 and rotor 2 are 994.88 rad/sec and axis of rotor 1 and 
rotor 2 are 994.88 rad/sec respectively, while their area moments of inertia are 
3.976x10"® m"* and 5.2x10"®, respectively. The naoduins of elasticity, E is 
2.07xl0"N/m^ With reference to Fig. 2.8, for roror 1, one gets 

{S}|=[f]2[i>],[F],[?liMo P-31) 

where, {S}^ is the .state vector to the of station 1. Continuing further, the 
following matrix relations can be obtained 

IS' ■ 



{S>f =[Ph{S)^+{R,) 

{S}‘ 

{S}'=tP],{S}f+{J!2) 

{S}’ =(?]7m6[i>]«[F],{S>* 

{S>? +{«2}] 

{S}? =[P],[F]6m<;[F]5[[i’]5[^l4[-f'l4m3{«>? +!^2>1 

{S)7 =m,[F]6m6[^']5[[^']5m4[f]4m3[m3(S}3‘ +W» + {«2}] (2-32) 

[P]i and [PJy are modified point matrix to include the effect of spring stiffness. 
Similarly, for rotor-2, one gets 

Ms =m,{s}s -«} 

Mfl = [i"]!! [2’], {S}* 

Mf2=[n2Mf2-M} 

(S}f2 = [i’ll2[-F]l,[niimio[2’1.0m»M<.(f'l8Mf -{*2} 

Mf2 = [2’].2(^']i.[f]u[2^],»[2’llot^'kI^'!9msfM.{S}a -<«.}1-{J!2> (2.33) 

The bearing reaction vectors are 
{«,} = {0 0 0 
{/!,> = {0 0 0 

Employing equations (2.34) and the following boundary conditions in equations 
(2.32) and (2.33), 

v.r. .I-.’/''' i' '■ 

A/,,, =jl/,,7 =2/,J =y,i2 =P , 1 ; . (2.35) 

ni=n7=n«=ni2=o 

the eigenvalue problem of equation (2.29) can be set up. 

The algorithm for solving the e^nv^^K prq^ip 1ms been written in MATLAB to 
compute the natural frequencies and mode shapes. The results obtained are given in 



Table 2.1 along with those obtained by Taipeng (1988) from the exact solution. The 
difference between the natural frequencies is less than 0.01%, which serves to 
validate the transfer matrix procedure and MATLAB algori thm 

Table 2.2 Natural frequencies of the dual rotor system (Case Study 1) 


Mode* 

Exact 

(Rad/Sec) 

Program 

(Rad/sec) 

B 

174.716 

174.72 

F 

400.436 

400.45 

B 

518.674 

518.68 

B 

632.988 

632.99 

F 

894.410 

894.42 

B 

1349.958 

1350.00 

B 

1397.403 

1397.40 

F 

1432.092 

1432.10 

B 

1805.751 

1805.80 

F 

1829.741 

1829.80 

B 

1868.237 

1868.20 

F 

21 19.298 

2119.30 

B 

2292.203 

2292.20 

B 

2637.412 

2637.40 

F 

3056.654 

3056.70 

F 

3422.186 

3422.20 

F 

3721.947 

3722.00 

F 

4752.938 

4752.90 


* B-Backward Whirl; F-Forward Whirl 


2.2 Finite Element Model 

Modeling of multi-spool rotors was carried out usii^ the general-purpose finite 
element software NISA v-7.0. This software is being used in a cost-effective manner 
in a number of rotordjmamic applications. However, its limitations, while modeling a 
dual rotor, with shafts rotating at difife^^t speeds, became evident during the course 
of this study. The steps to be followed while modeling a rotor in NISA are briefly 
described here. . ,, . 

NISA employs DISPLAY-HI as a j^ and post processor, and NISA-ADVANCED 
DYNAMICS to csdculate and vedors of rotor system. Modeling 



process in NISA can be divided in two broad steps - Geometric Modeling and Finite 
Element Modeling. 

2.2.1 Geometric modeling 

The geometry of the dual rotor, in the present study, has been represented by 
employing a solid primitive cylinder and hyperpatch. For geometric modeling macro 
provision of NISA is used. It is a set of statements (program) which can be executed 
in DISPLAY III to model geometry in terms of variables and thereby can be used to 
create parametric models by inserting different values for the variable. 

2.2.2 Finite Element modeling 

3D solid hexahedral elements (Fig. 2.9) - NKTP = 4, DOF = 3 (UX, UY, UZ) are 
used to model the dual rotor . The discs are also modeled using 3D solid hexahedral 
elements. 3D translational spring elements (NKTP=17,DOF=3 UX, UY, UZ) are used 
to model intershaft and end bearings. All required geometric and material properties 
are defined through the NISA II ‘forms’ structure. 


1 2 



Figure 2.9 Hexah^ral dement with Nodes 1 to 8 
2.2.3 Case Study 2 

The modeling in NISA and computation of natural fiequencies and mode shapes is 
illustrated through the example of the dual rotor chosen by Rajan, Nelson and Chen 
(1986). The rotor has also been an^yzed through the transfer matrix method 
described in earlier sections, for verificatkffl. of mode shapes. The rotor is shown in 
Figs. 2.10 (a), (b) and the physical and gemndric data are given in Table 2.3 (a)-(c). 


II 



Table 2.3 (a) Rotor Subelement Data (Case Study 2) 


Shaft 

No. 

Element 

No. 

Subelement 

No. 

Length 

(10'^) m 

Inner 

Radius 

(10'^) m 

Outer 

Radius 

(10'^) m 

1 

1 

1 

7.62 

0. 

1.524 


2 

1 

8.89 

0. 

1.524 


2 

2 

8.89 

0. 

1.524 


3 

1 

7.62 

0. 

1.524 


3 

2 

762 

0. 

1.524 


4 

' 'J s, 

5.08 

0. 

1.524 


5 

1 

5.08 

0. 

1.524 

2 

)■ 

7 

1 

5.08 

1.905 

2.54 


g 

C’-' ”1 

7.62 

1.905 

2.54 


8 


M2 

1.905 

2.54 


,:9 

a.:, 

_ . IT — 

5.08 

1.905 

2.54 


i E = 206.^10^ i \ * - P = 8304 Kg/m^ 


Table 2.3 (b) Rotor Concentrated Disk Data (Case Study 2) 


Shaft 

Station 

Mass 

Polar Inertia 

Diametral Inertia 

No. 

No. 

(Kg) 

(kg-cm^) 

(kg-cm^) 

1 

2 

4.904 

271.2 

135.6 


5 

4.203 

203.4 

101.7 

2 

8 

3.327 

146.9 

073.4 


9 

2.277 

097.2 

048.6 

Table2.3 (c) 

Rotor Bearing Data (Case Study 2) 

Shaft 

Station (a) 

Station (b) 

Stiffness 

No. 




(N/cm) 

1 


1 

0 

262,795 



4 

10 

087,598 



6 

0 

175,197 

2 


7 

0 

175,197 


(0 denotes a connection to the rigid base) 


The NISA model of the rotor, along with the co-ordinate system used are shown in 
Figs. 2.10 (a), (b). It needs to be mentioned that simulation of the shaft-disk assembly 
needs the additional procedure of simulating the interference fit between them. This is 
done by making the nodes, on the inner periphery the hole of the disc, lie on their 
counterparts of the shaft’s outer periphery. The nodes are then merged with specified 
tolerance. The tolerance decides the type of fit, which is required to be obtained in the 
assembly. Close nodal merging tolerance of 0.001mm were ensured during assembly 
in order to simulate the interference fit condition of real life rotor system. A major 
limitation posed by NISA software is that it does not provide for rotation of 
components. Effects like gyroscopic couples can be simulated only through specific 
application of torque, which needs to be computed sepmately, for every component 
and rotating speed. During the present study, free vibration solution in NISA, does 
not incorporate gyroscopic influences. The NISA solution however provides good 
estimates of the natural fiequencies, since gyroscopic couples are hi^er order effects. 
The mode shape aHJroximation is also good and provides with ftiree-dimensional 
facilities for motion visualization. The natural frequencies of the rotor obtained 
through this FEM model along with those from Rajan, Nelson ai^ Chen (1986) and 

the transfer matrix method are given in Table 2.4. 
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Figure 2.11(a) 


:e system for 3D FE model 
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T able 2.4 Natural Frequencies (rad/s.) of Dual Rotor (Case Study 2) 


Rajan et al 

(1986) 

T M Method 

FEM (NISA) 
(Solid Element) 

(FEM) NISA 
(Beam Element) 

755.0 

771.0 

755.9 

749.42 

1531.0 

1593.0 

1445.5 

1540.92 

2244.0 

2220.0 

1949.5 

2261.40 


It can be observed that the beam element model in FEM gives natural frequency 
values, which are closer to the Transfer Matrix model of the present study and that of 
Rajan et al, since the beam element simulates similar degrees of freedom as those in 
the Transfer Matrix models. Solid element model however provides the advantage of 
visualization of stress distribution, resultant displacements and contour plots. The first 
three mode shapes of the rotor corresponding to transverse vibration are given in Figs. 
2.12 and 2.13. Figs. 2.12 correspond to the beam element model, while Figs.2.13 
pertain to the solid element model. The mode shapes from the two models can be seen 
to be similar. The degree of flexure in the outer rotor in all three cases appears to 
significantly less than that of the inner rotor on which it is supported through the 
intershaft bearings. In the second mode, the nodal point for the inner rotor is close 
intershaft bearing number 2 (on the right), while in mode 3 the nodal point has shifted 
towards intershaft bearing number 1 . This is also reflected in the contour plots, which 
show maximum stress at the junction points of the two rotors. 


Tht critical speeds of this rotor, for forward and backward synchronous whirl with 
inner and outer rotor, and the corresponding mode shapes have been obtained by 
using transfer matrix method. These results are listed in Table 2.5 (a), (b). 


Table2.5 (a) Critical Speeds (rad/s) - Synchronous whirl with inner rotor 

(Case Study 2) 

F orward Whirl Backward wMrl 

Raian et al Present Study Rajan etal 
(1986) T.M.M. 

863 895 660 676 

1606 1692 1423 1458 

2283 2269 2125 2090 _ 
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Figure 2.12 Modeshapes (case study 2) - beam element model 


26 
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Figure 2.13(a) I Mode (Case Study 2)-SoKd Element 

(I)ModeShape 
( n ) Contour Plot 


27 


HI - OEOhCTRY HOUELING SYSTEM <8*0, 0> PRE/POST MCttMJLE 


tKMt SHflPC PLOT 
m »SF» 1.04£-€2 

Hom HO.* uses 
sa&E » i.e 
mppm scaiMO 




Figure 2.13(b) n Mode (Case Study 2)-Solid Element 

(I) Mode Shape 

( II ) Contour Plot 
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Figure 2.13(c) m Mode (Case Study 2)-Solid Element 
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FW- Forward Whirl BW = Backward Whirl 


Figure 2.14(a) Mode shapes corresponding to critical speeds (Case Study 2) 

( Synchronous whirl with inner rotor ) 


30 












0.2 0.3 0.4 

Length of Rotor 





FW = Forward Whirl 


BW = Backward Whirl 


Figure 2.14(b) 


Mode shapes correspouding to critical speeds (Case Study 2) 

( Synchronous whirt wiA outer rotor ) 
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Table2.5 (b) Critical Speeds (rad/s) - Synchronous whirl with outer rotor 


Forward Whirl 

Backward Whirl 

Rajan et al 

Present Study 

(T.M.M.) 

Rajan et al 

Present Study 
(T.M.M.) 

882 

851 

687 

704 

1584 

1664 

1462 

1505 

2274 

2259 

2175 

2140 


There is close resemblance between the results from the component synthesis method 
of Raj an et al and those from the transfer matrix method of the present study. The 
error can be reduced by increasing the number of stations in the transfer matrix 
model. The mode shapes corresponding to the critical speeds are shown in Figs. 2.14 
(a), (b). It can be seen from these shape acquired by the rotor during forward and 
backward whirl is similar in a particular mode. However, the relative separation 
between the two rotors are different and is higher in backward whirl for the case 
investigated. The trend is similar for synchronous whirl with outer rotor and inner 
rotor. Further studies are required to describe the phenomenon in more detail. 

Critical speeds have also been obtained for non-synchronous whirl of the rotor. These 
results are shown in Fig. 2.15 for the speed ratio of 1.5. The synchronous whirl 
critical speeds can be obtained from the intersection of straight lines AB and AC with 
the curves, for inner and outer rotors respectively. 


wout •Rotational speed of outer rotoi' (rad/s) 



Figure 2.15 Natural frequencies (NF) and critical speeds of dual rotor 
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CHAPTER 3 


UNBALANCE RESPONSE OF DUAL ROTOR 

Residual unbalance is one of the most significant causes of excessive vibrations in 
rotating shafts. Unbalance primarily arises due to material inhomogeneities, 
manufacturing defects, keyways, slots etc. In addition deterioration due to wear and 
tear, thermal bending, process dirt collection etc. during operation also give rise to 
eccentricity between geometric and gravity centres, causing unbalanced forces. These 
can be removed by proper balancing procedure. It is important to determine response 
of the rotor due to unbalance, in order to study its dynamic behavior and set limits for 
safe operation. 

3.1 Transfer Matrix Procedure 

The extension of the transfer matrix procedure, discussed in the previous chapter, for 
unbalance response analysis of multi-spool rotors is described here (Rao, 1 996). 


3.1.1 Field Matrix 

With reference to Fig 2. 1-2.3 and accounting for bending in both x-z and x-y planes, 
the state vector {S} can be expressed as 



where 



V 


-w 




e 

My 


-Vy 

1 y j 


.K. 


It can be readily shown that 

where [F] is the same as defined previously by equation (2.8). Relations (2.8) and 

(3.3) can then be combined to give the overall smte vector relation as foll^^^ 
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m [Ollte)!* 
_[0] [f]Jl{s,}J.__ 


(3.4) 


The unbalance mass in>'-z plane is shown in Fig 3.1. The eccentricity, e, at station i 
makes an angle /? i in yz coordinates fixed to the rotor. The unbalance then is 


u 


y,‘ 


m-fi, cos/?,.; 


= »J,e,.siny0,. 


(3.5) 



Figure 3.1 The ith unbalance mass 


In general the response in the y-direction will be different from that in the z-direction. 
Further for anisotropic bearings the response in both the x-y and x-z planes will 
consist of sine and cosine components. Therefore splitting {Sz} and {Sy} further as 



and adding the identity 1 s 1, to the 16 equations of total state vector for the purpose 


of response calculations, the overall field equation becomes 


McY 

L 

{S.s} 





I 


[F] 

£0] 

[0] 

[0] 

{0} 

£0] 

IF] 

[0] 

[0] 

{0} 

[0] 

[0] 

[F] 

[0] 

{0} 

[0] 

[0] 

[0] 

[F] 

{0} 

{0}’^ 

(0}^ 

{0}" 

{Of 

1 





{SJ 

< 


Ji 



which is 


(3.9) 
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{S}f =[F],{S}f_, 

where {S} and [i^] are modified state vector and field matrix respectively. 


3.1.2 Point Matrix 

The fi'ee body diagram for the ith unbalance mass in the x-z plane is given in Fig 2.3. 
Noting that and similar expression can be written for jc-y plane. 
yR _ j/i _ fn.co\ - co^u , cos cot + Q}^u^,sin(jot 

'/ 2 /’ /’ (3-11) 

Vy i = Vy t - mfi) w,. - (O Uy . co^cot + co u^ . sin cot 
the relationships across the fth mass are written as 

\[P] [0] [0] [0] Kjife}]" 

[0] [P] [0] [0] {m,,} {5,,} 

^{Sya}\ = [0] [0] [P] [0] {S',,} (3.14) 

{5,,} [0] [0] [0] [P] {m,,} {S',,} 

1 . {0}'’ {0}^ {0}^ {0}" 1 1 

K J I L J i 


where [P] is the point matrix defined earlier in equation (2.10) with critical speed p 
replaced by rotational speed o , and 



Equation (3.14) can be written as 


{S}f =[?],{«>? 


3.1.3 Bearing Matrix 

The stil&ess and Hamping properties of bearings are incorporated in the bearing 
matrix. Refering to the equilibrium relations across a bearing (Fig 3.3) the bearing 
transfer matrix is 

Wf =[«],{«>? 
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where the non-zero elements of [S] are 


5,,,=l,r = 5 = l,2 17 

^4.1 =^8,5 


■^ 12,9 -^ 16,13 ~ ^yy -^ 4,9 “ -^ 8,13 ~ ^zy 

^ 4.5 ~ ~^ 8,1 ~ ~^Z2^ “ ~^S,9 ~ ^zy^ 


•®4,13 ~ -^8,9 ^ zy ^ 

-®12,1 ~ ~-®16,5 “ ^yz 


■®12.5 " -^16.1 ~ ^yz 

^2,13 “ “^6,9 ~ ~^yy^ 


(3.16a) 


The overall transfer matrix of the rotor system mounted on fluid film bearings at both 
ends 


<S>^,=[SU[C/)[S]o{S>5 (3-17) 

{S)«l=[^»]{S>J (3.18) 


Where \Um\ is the overall transfer matrix of size 17x17, which includes bearings 
matrix also. For rotor with bearings at its both ends (the shear forces and bending 
moments are zero at both ends), one obtain the following 

equation: 


^3,1 “3,2 “3,5 “3,6 “3,9 “3.10 “3,13 “3,14 



“3,17 

“4.1 “4,2 — • 

0 c 


“4,17 

“7,1 

- W , 


“7,17 

“8,1 

0 s 

y — - ^ 

“8,17 

> 

“H,l 



“11,17 

“12,1 



“12,17 

“l5,l 

^s 


“15.17 

_“.6.1 

. . 

0 

“16,17. 


The above system of linear equation can be solved to determine w^,6cetc., at the 

zeroth station, which will define {S}J. The state vectors at all other stations are then 
obtained appropriate transfer matrices to get the unbalance response of rotor. 
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3.1.4 Modifications for Dual Rotor 


In a dual rotor presence of intershafl bearings requires modification in the state vector 
with field and point matrix, since the dynamic response of each rotor not only 
depends upon its own imbalance but also the dynamic behavior of other. Therefore, 
each quantity of the state vector requires to have two components, one corresponding 
to inner rotor speed and other to outer rotor speed. If the suffixes m and n represents 
the inner and outer rotors respectively, - inner rotor speed, O)^ - outer rotor speed 
etc. Then for x-z plane we have 


'W = 'Wcfyi coscojfit + sinui/^jt cosui^t sin<U;^r 

^ +^5771 sin cos sin 19) 

My = Mycm coscOfn^+Mysm sincOfn^+Mycn cosoint+^ysn sinui^r 
Vz ~ ^zcm C'OSOmt'^^zsm cosui^f cosui^f 

Similar expressions can be written for the x-y plane and the state vector is modified to 
contain 33 quantities 


{S} = K,> {S,,} is,.) {S,} if 

where 


{^zm } 






^cm 


^cn 

d^m 


^cn 


<l>cm 


<t>cn 

Mycm 


Mycn 


Mzcm 


^zcn 

'XT 

y 

^ zcm 


V 

^ zm 


^ycm 


^ yen 

r 

V ^ 

-'^sm 


^sm 


^sn 

^sm 


^sn 


^sm 


^sn 

^ysm 


Mysn 


Mzsm 


Mzsn 

V 

\j zsm . 


V 

y zsn . 


V 

^ ysm 


J^ysn 


(3.20) 


(3.21) 


Now it can be readily shown that corresponding field matrix is given by 
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[^■1 


33x33 “ 


0 

0 

0 

0 

0 

0 

0 

0 


0000000 0 ' 
[■^1 0 0 0 0 0 0 0 

0 [F] 0 0 0 0 0 0 

0 0 [F] 0 0 0 0 0 

0 0 0 [F] 0 0 0 0 

0 0 0 0 [F] 0 0 0 

0 0 0 0 0 [F] 0 0 

0 0 0 0 0 0 [F] 0 

0 0 0 0 0 0 0 1 


(3.22) 


where [F] is the field matrix derived in section 2.1. Also the point matrix can be seen 
to be readily modified to 



[p] 


33x33 “ 


0 

0 

0 

0 


0 

[Pn] 

0 

0 

0 


0 0 {mji 

0 0 

{PJ 0 {niyj 

0 [Pn\ {^yn) 

0 0 1 


(3.23) 


where 



0 






■[^] 

0 


0 




p-^fl 


For inner rotor {m ^ } and } will be null matrix and {m^},{mym } will be given by 
equation (3.24). 

{m™} = {0 0 0 uyj 

{m^} = {0 0 0 uyaOduyf (3.24) 

and [Pj is given by equation (2.10). If the gyroscopic couple is to be taken in 
account, [F] in the above equation can be taken from the equation (2.24). Similar 
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expressions can be written for the outer rotor. Also, the bearing matrix gets modified 
to 


[Be^i] 

0 

\.B^mz2\ 

0 

O' 


0 

[Be„,i] 

0 

{Be„,2\ 

0 


[Be„y2] 


[Be^yx] 

0 

0 

(3.25) 

0 

l-B^nyll 

0 

[Be^^] 

0 


0 

0 

0 

0 

1 



where non zero terms of sub matrices are 

= = U 8 

^ ^zz^m 

[Be^2]m = [38^2^5) = -K^ 

[5em2](4,5) = C^o)„ 

[Be^2W,l) = -C^a2^ 

Matrices [Bcmyi] and [Bemyi] are obtained from the above by replacing subscripts z 
with >' and y' with z. Similarly one can get the other sub matrices used in bearing 

matrix. 


3.1.5 Junction conditions 

With reference to Fig. 3.3 if [A], [B], [C] and [D] represent the overall transfer 
matrix considering all the stations between the starting station zero, for the A, B, C 
and D shaft sections and the junction point where the intershaft bearing is located. 


then 

-[^]{ 5 }o^ {S}oB =[B]{B}mj 
-[C]{5}oc =^[B]{S}^j [BeHBeHn 
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Figure 3.3 Junction Conditions 


At the junction of the inner rotor, one gets 

Similarly at the junction of outer rotor 

{S)%=[IWf^jHBe}{S)‘-„j-\.Bemij (3 28) 

Substituting (S}Bj and {S} * 6om the above in the second and fourth equation of 
(3.26), 

CToi = p 29 ) 

(S>0D =P][[f]+tSe)]{S}J; -[II][Se]{S}^j 

With the help of first and third equations of (3 .26), the above can be written as 

(330) 

Premultiplying the first equation by [B]~^ and ^ond equation by [D] ^ the above 

equation can be written in one matrix equation. 
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(3.31) 


[AJ [BJ [CJ [0]' 
lA] [0] [CJ [D„] 


66j:132 


{5}o. 

{^}0B 

{■5}oc 

{5}oz) 


= 0 


132x1 


where 


[AJ = [[I] + [Be]][A] 

[BJ = -[Br^ 

[C„] = -[Be][C] 

[A„] = -[Be][A] 

[CJ = [[/] + [5J][C] 

[D„] = -[D]-‘ 

If the intershafl bearing matrix is unit matrix, it can be seen that equation (3.31) 
reduces to two uncoupled equation for the outer and inner rotors. For free end 
boundary conditions at the four ends of dual rotor system, all moments and shear 
forces are assumed to be zero, i.e., at these ends, 

^ ym ~ ^ yn ~ ^ zm ~ ^zn 

V =F =V =V =0 
'^ym ^ yn zm zn ^ 

The unknown quantities are M>,d,v,^ at the four starting stations. For each starting 
station there are 16 unknowns 

T 

•••••A^sn^^sn^QA 

and similarly {S}ot,{5}oc and {S}m unknown quantities at the starting 

stations of all four section. Since My,V^ etc, terns in the vector of equation are zero. 

Columns 3,4,7,8,11,12, ...31,32; 36,37...64,65; 69,7€...97,98; 102,103.....130,131; 
can be dropped in the 66x132 matrix of equation. The 33“^^, 66 , 99 and 132 
elements of vector are unity and hence whai multiplied with the corresponding 
column quantity in the matrix, give rise to a constant quantity independent of any 
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state vector dynamic quantity. These four terms can be added and taken to the right 
hand side of the corresponding equation. 33^^* and 66^^ rows of the equation (3.31) can 
be deleted as they simply give identity equations. Then the equations (3.31) reduce to 

[^]64x64{So}64xl = {^)mx\ (3.32) 

where 


{5o} = 


'{S}0A 

{§}oc 

{S}0D 


^ ‘ ^mi,33 ^mi,33 ^mi,33 


(3.33) 


and A(i,j) are the elements picked up from the matrix (3.31) as described in the 
above manner, by dropping columns 3, 4,7, 8 etc., and rows 33 and 66. 


The linear system of equation (3.32) can be solved to obtain {Sq} in (3.33). From 
there, the state vector quantities at any station can e determined using the appropriate 
transfer matrices. 

The response of a straight rotor is characterized by siqjerposition of two perpendicular 
harmonic motions with the same angular frequency, which gives rise to an elliptical 
motion. However, for dual rotors, the unbalance response comprises of two different 
frequency components, belonging to the inner and outer rotors running speeds. In 
such a case, the sum of the different harmonic functions will be a periodic function if 
the ratio of the two frequencies is a rational number, and beating phenomenon will be 
observed with the an angular frequency given by die largest common divisor of both 
frequencies. (Rossing et.al. 1995), The period of such a beating phenomenon is 
T = 2n:{p-q)l{a>n-a)rft), (3.34) 

where p at yl q are the smallest integers suchtiiat p! q — <o„l The superposition 
of the response components wand v leads to a motion characterized by a closed 
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periodic orbit, with a number of loops. The number of such loops is given 
n = {Max[(On,0}yn]}{p-q)l{o)n-o)^). (3.35) 

The algorithm for unbalance response computation has been developed in MATLAB 
and verified with the results of a case study taken up by Gupta et al (1988). 

3.2 Case Study 3 

The extended transfer matrix method for dual rotor described above, is illustrated 
through the idealized dual rotor model mounted on isotropic bearings as shown in 
Fig.-3.4. The shaft diameter is 1.5 cm and the Young’s modulus is taken as 210x10^ 
N/m^. The unbalance on each mass is in phase with the other and is equal to 0.001-kg 
m. The speed of inner rotor is 1200 rpm and outer rotor is 1000 rpm. As discussed 
earlier, the response of each rotor, will be affected by the unbalance of its own as well 
as the other rotor. The total response at any speed to be the sum of two harmonic 
components, in present case with 1200 rpm and 1000 rpm frequencies, giving rise to a 
beating phenomenon. The whirl orbit for disc land 2 of both rotor at specified speed 
ratio is shown in Fig.3.5. 



Figure 3.4 Idealized dual rotor model 

The first two natiral fiequencies of flie present dual rotor model can be computed to 
be 238 tad/s (2273 q.m) and 539 rad/s (5147 rpm). The unbalance response has been 
computed for two sets of inner and outer rotor speeds. The first set is chosen as 2500 
rpm for the inner rotor and 2000 rpm for the outer rotor. This set lies close to the first 
natural frequency. The other set of inner and outer rotor speeds is chosen as 5500 rpm 
and 5000 rpm, so that the operation is closer to the second natural ftequency. For the 
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Figure 3.5 Whirl orbits at speed ratio of 1200 rpm/1000 rpm 
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Whiri orbN: at disc 1 of rotor 1 



X 10'^ Whirl orbit at disc 2 of rotor 1 



WNrI orbit at disc 1 of rotor 2 



Whirl orbit at (fisc 2 of rotor 2 



I 


Figure 3«7 Whirl orbits at speed ratio of 5500 rpm/5000 rpm for Dual rotor 


i' 
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first set of speeds p/, = 5/4, r = 0.12 s and n = 5. For the second set. we have 
p! q 11/10, T 0.12 s and n — Il.The unbalance response in the two cases is 
shown in Figs. 3.6 and Figs 3.7, respectively. It can be seen that the orbits traced out 
by the inner and outer rotors at the two speeds are different in many ways. This 
dynamic behaviour would be dependent on a variety of design factors and operating 
speed regimes. An attempt has been made in the present study to investigate one such 
design parameter, i.e. the cross-coupling effects of the intershaft bearing. T his is 
discussed in the next section. 

3.3 Influence of Cross-Coupling Stiffness of the Intershaft bearing 

The influence of the value of the cross-coupling stiffness of the intershaft bearing on 
the unbalance whirl orbit has been investigated by defining the cross coupling 
parameters 

M\ ^yz ^ ^yy 

M2 ~ ^zy ^ ^zz (3.36) 

For the present study the two parameters have been taken to be equal, i.e. 

= A2 = A • Retaining all parameters, from the previous section, including 
kyy = =1.0 MN/m , the values of cross-coupling terms ky^ are varied to 

obtain whirl orbits for p = 0.0, 0.5, 1.0 . These orbits are shown in Figs. 3.8 and 3.9 
for two speed ratios used earlier, namely 2500rpm / 2000rpm and 5500rpm/5000rpm. 
It can be seen from these figures that cross-coupling values play a more critical role 
in determining the unbalance locus of the rotor, for the second set of speeds. The first 
set of speeds is in the vicinity of the first mode of the rotor, whereby both shafts are 
moving in phase. The second set of speeds is in the vicinity of the second mode, 
whereby the two shafts move in op^wsite phase and more transfer of energy takes 
place through the intershaft spring. Role of die cross-coupling parameter p is mote 
critical here as evidenced in Fig.3.9. Also, for both sets of speeds increase in cross- 
coupling results, in increase of the skewness of the orbit. A reduction in the orbit 
amplitude is also visible, with increase in p - a result of increased system stiffness. 
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No attempt has been made in the present study, to investigate the influence of 
intershaft damping. A comprehensive analysis of the dynamics of multi-spool rotors 
and its dependence on various shaft geometry parameters, speeds and speed ratios, 
bearing locations, stiffness and damping parameters requires more exhaustive 
simulation. 
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CHAPTER 4 


DYNAMIC ANALYSIS OF AEROENGINE ROTOR 

In this chapter, some aspects of the dynamics of an aeroengine rotor are attempted to 
be analysed. A typical aeroengine rotor with two shafts is chosen for study. The inner 
shaft carrying the LP stages of the the compressor and turbine at its two ends and 
supports the outer rotor in its central portion, through two intershaft bearings. Free 
vibration analysis is carried out using - both the transfer matrix procedure and FEM 
package. Residual unbalance response is computed using transfer matrix procedure. 

4.1 Stage Data 

The aeroengine rotor is schematically shown in Fig. 4.1. The outer rotor carries the 
HP stages of the compressor and turbine. The compressor comprises of 5 LP stages 
and 6 HP stages, while the turbine comprises of one LP and HP stage, each. The 
various stages are modeled as discs. Most of the data is mentioned on the Fig.4.1 
itself. Some other relevant features are listed in Table 4.1 . 


Table 4.1 Aeroengine rotor Disc Diameters 


Component 

Name 

Stage No. 

Outer Diameter 

(mm) 

Inlet vane 


800 

L.P. Compressor 

I 

800 


2 

720 


3 

640 


4 

560 


5 

480 

H.P. Compressor 

1-6 

400 

H.P. Turbine 

1 

800 

L.P. Turbine 

1 

^ •! /vi) 'KT/— 

575 

^2 iTrrU 


Young’s modules (E) 210x10^ N/m^ Density 8304 Kg/m 


(material properties of all parts are taken as same in die Resent model) 
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Figure 4.1 AEROENGINE 
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Figure 4.2 3D Finite Element Model of Aeroengine Rotor 
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4.2 Bearing Data and Stiffness Computation 

The bearings supporting the inner and outer rotors and the intershaft bearing, in the 
aeroengine rotor are all rolling element t3^es. Data for the typical bearings chosen for 
this case study are given here (refer Fig. 4.1). 


The bearings have modeled as two translational springs, one in each y and z direction. 
Cross-coupling effects and damping are neglected, since these effects are small in 
case of ball bearings. Stiffness computation is done by employing the procedures 
outlined by Harris (1966), based on the extension of Hertzian contact theory. The 
procedure is briefly described here. 

The total elastic force at the point of contact of the Mi ball with the inner or outer 
races is expressed as 

h = + +Tsin77jt)^^^ (4.1) 

and its projection along the line of action of applied force is 
F; ==i^„(g + xcos7*+ysin;7*)^^^cos77^ (4.2) 

where g is the radial preload or the radial pre-clearance between the ball and the races 
and X and y are the displacements of the moving ring in the direction of radial load 
and perpendicular to the direction of radial load respectively, is the angle between 
the line of action of the radial load (direction of displacement of the moving ring) and 
the radius passing through the center of the kth. ball. K„ is a coefficient of 


Table 4.2 Bearing Specifications 


Bearing 

No. 

Type of 
tearing 

No. of 
rollers/balls 

Roller/Ball 

dia(mm) 

Inner 
race dia 
(mm) 

Outer 
race dia. 
(mm) 

Roller 

Length 

(mm) 

Radial Clear - 
- ance (mm)xlO' 

Ninntrl 

N inner 2 

Nonte*! 

Nouter2 

Nintarshl 

N f 

Roller 

Roller 

Ball 

Roller 

Bdl 

Roller 

22 

30 

29 

30 

21 

36 

09.5 

11.0 

35.0 

15.0 

24.0 

06.0 

080.0 

120.0 

225.0 

150.0 

120.0 
080.0 

140.0 

215.0 

405.0 

270.0 

215.0 

140.0 

25 0 

40 0 

55.0 

25.0 

0.085 -0.115 
0.095-0.130 
0.150-0.190 
0.030 - 0.080 
0.110-0.160 
0.045-0.142 
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proportionality depending on the geometric and material properties of the bearing. 
The coefficient of proportionality is estimated by the method suggested by Harris 
(1966), as follows. Writing equation (4.1) as in general for ball bearing. 

^ ~ exponential value is 1.1 1 for roller bearings) (4.3) 

where S is defined as a deflection. 

and since he total normal approach between two raceways under load separated by a 
rolling element is the sum of the approaches between the rolling element and each 
raceway, i.e. <5^ = Sj^. + 

Therefore 

K,=[lK(XIK,f"HVKj'’)r (4.3) 


F = 

(4.4) 

For steel ball-steel raceway contact. 


=2.15xlO®Sa-“‘^^(^*)-3^2 

(4.5a) 

and for steel roller and raceway contact. 


=0.786x10^/®^’ 

(4.5b) 


where Ecx is curvature sum and 5* is dimensionless deflection or contact 
deformation and is a function of the curvature difference da and / is roller length. 
(Hairis, 1966). Curvature sum for ball-inner ring raceway and for ball-outer ring 
raceway can be calculated from the following set of equations. 

1 


' f, 1-r 


10 - 0 =— ( 4 - — --^) 
D /o 1 + r 


£)cosa 


da. 




1+iL 

fi i-r 

. 1 2y 

4 — .-4.. 


/ = - = - 
2(p D 

1 2y 


da„ 


/o 1 + r 

2y 


fi 1-r 


4_ J-. 


(4.6) 


1+r 


TTte total elastic force in the direction of the applied force is (refer equation 4.2) 
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(4.9) 


F=iFk 

k=l 

where n is the total number of balls in the bearings. 

Using the condition of zero elastic force in the direction perpendicular to the elastic 
load, the deformation y, perpendicular to the radial force line is expressed as 

i [g + xcosTJkf^^sinTji^ 

y = — — sin^ (4.10) 

2:[g + xcos7j' ^ 

The bearing stiffness is determined as the function of the deformation x as 
dF 

k{x) = ~ (4.11) 

ox 

Stiffness values thus calculated, for various bearings on the aeroengine rotor are listed 
in Table 4,3. The bearing is assumed to be isotropic, with identical stiffness in all 
radial directions. Cross-coupling influences in rolling element bearings are generally 
negligible. Damping is also small in these bearings and is ignored in the present 
model. 


Table 4.3 Bearing Stiffness 


Bearing Number 

Bearing Stiffness 

(N/mm) 

Ninner 1 

2.6240x10'’ 

N inner 2. 

5.7600x10® 

Noutcr i 

0.1096x10® 

Nouter 2 

7.0480x10® 

Nintersh 

0.0915x10® 

Nintersh 2 

2.6240x10® 


4.3 Natural frequencies and Mode shapes 

The FE model of aeroengioe rotor, using 3D solid hexahedral elements is shown in 
Fig.4.2. Free vibration analysis is carried out using NISA software as well as the 
transfer matrix algorithm in MATLAB. A comparison between the two sets of results 
is given in Table 4.4. The natural ftequencies from the two procedures are fairly 
close, the difference primarily arising from the disoetisation process m Transfer 

matrix method. 
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Figure 4.3 I Mode of Aeroengine Rotor 


(a) Modeshape 

(b) Displacement Contours 

(c) Stress Contours 
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Table 4.4 Natural Frequencies (rad/s.) of Aeroengine rotor 


Transfer Matrix Method 

FEM model 

226.00 

226.76 

424.00 

438.06 

443.00 

489.58 

717.00 

572.16 


Fig. 4.3 shows the modeshape corresponding to the first natural frequency of 
transverse vibration, This figure also includes the contour plot of resultant 
displacement and stress in aeroengine rotor for &st natural frequency. Similarly the 
Fig. 4.4 and Fig 4.5 shows the modeshape corresponding to second and third 
transverse natural frequency respectively. Critical speeds of aeroengine rotor for 
forward and backward synchronous with inner and outer rotor are given in Table 4.5. 


Table 4.5 Critical speeds (rad/s) of Aeroengine rotor 


Synchronous Whirl with inner rotor 

Synchronous whirl with outer rotor 

Forward Whirl 

Backward Whirl 

Forward whirl 

Backward Whirl 

312 

172 

282 

186 

764 

295 

667 

325 

1104 

337 

781 

358 

1544 

563 

1255 

604 


4.4 Residual Unbalance Response 

The unbalance responses of aeroengine rotor at different speed ratios are obtained by 
extended transfer matrix method. Rotors in an aeroengine rotor are balanced at 
regular maintenance intervals. However, some residual unbalance still remains, 
depending on the quality of balancing achieved. Unbalance response has been 
computed here, for such typical residual unbalance, which is listed in Table 4.6. 
Whirl orbit for inner and outer rotors have been obtained for three different speed 
ratios. Fig 4.6 shows the whirl orbits at the speed ratio of 2500rpma000ipm. 
3500rpm/3000rpm and 4500 tpin/ 4000 rpm. The time period of whirl was found to be 
0.12 s, with 5,7 and 9 loops for first, second and third set of speed ratio respectively 
fiom equations 3.34 and 3.35. Hie whirl oibHs have been plotted at the fourth stage of 
LP compressor for inner rotor and at the third stage of HP compressor for outer rotor. 
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Table 4.6 Unbalance on Aeroengine rotor 


Component Name 

Stage No. 

Unbalance (miCi) 

Kg.m 

Angle of eccentricity 

(A) 

Inlet vane 

- 

0.0372 

20.0 

L.P. Compressor 

1 

0.0180 

10.0 


2 

0.0110 

5.50 


3 

0.0080 

4.00 


4 

0.0060 

3.00 


5 

0.0110 

5.50 

H.P. Compressor 

1 

1.200 

6.0 


2 

1.200 

6.0 


3 

1.200 

6.0 


4 

1.200 

6.0 


5 

1.200 

6.0 


6 

1.600 

8.0 

L.P. Turbine 

1 

0.063 

31.5 

H.P. Turbine 

1 

0.880 

4.50 


The whirl orbits for the various stages along the rotor length at the speed ratio of 
2500rpm/2000rpm are shown in Fig. 4.7 (a)-(b). The plot however, does not correctly 
represent the relative orientation of the orbits of the various stages. The dynamic 
displacements at various stages, along the rotor, at one particular time instant, have 
been shovm joined by straight lines, in the figure. From Fig 4.7 (a)-(b) one can 
visualize that at given speed ratio of 2500rpm/2000rpm, which is close to first natural 
frequency of aeroengine rotor the mode of vibration of inner and outer rotor is very 
close to first modeshape of aeroengine (Fig 4.4) .Better post-processing of results is 
required for a more realistic spatial representation of the rotor response. The exercise 
is however, sufficient, to highlight the fact that aeroengine rotor dynamics, is a 
function of a variety of design and operating parameters and a rigorous parametric 
study is essential to draw meaningful inferences. 
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CHAPTER 5 


CONCLUSIONS AND SCOPE FOR FUTURE 

WORK 


Models and algorithms have been developed during the course of this study for 
dynamic analysis of dual and multi-spool rotors. The basic features of such rotors are 
rotating shafts, which support other rotating shafts through intershaft bearings. This 
arrangement, while improving upon the compactness of the installation makes its 
dynamic behaviour more complex by introducing additional coupling influences of 
intershaft bearings. The transfer matrix models account for such inter-component 
interactions in a simple and efficient manner. However, the accuracy of results is a 
function to the discretisation adopted. Modeling of rotating components is difficult in 
most finite element software packages. The difficulty is compounded in the case of 
dual rotors, with shafts rotating at different speeds. Keeping these limitations in view, 
modeling was attempted in NISA software package in order to appropriately model 
intershaft comiections and obtain better visualisation of the complex mode shapes. 


These models could be successfully developed. The results for standard problems 
available in literature were obtained to validate the models. Gyroscopic effects were 
included in the transfer matrix algorithms and critical speeds have been calculated for 
forward and backward whirls synchronous with speeds of different rotors. A chart for 
nonsynchronous whirl as function of rotating speeds is also generated. Modeling and 
analysis has been carried out also for an aeroengine rotor. Stiffiiess calculation for 

typical rolling element bearings has been carried out. 


The present models establish the framework for detailed parametric studies on dual 
rotor dynanucs. Influence of various intershaft stiffiiess parameters needs to be 
investigated in more detail. This can include both direct stiffness and eross-stiffiiess 
terms. Influence of damping can also be smdied wift the help of pre^t computer 
programs. Influence and location of other bearings on the mner and outer shafts can 



be studied m terms of nondimensional parameters. Rigorous studies can also be 
conducted to investigate the coupling influences at different speed ratios. 
Experimental investigations can also be carried out, by building simple laboratory 
models, to validate the computational results and investigate other phenomenon like 
instability in dual rotors. 
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